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We show that the two-dimensional, nonlinear Schrodinger lattice with a saturable nonlinearity 
admits periodic and pulse-like exact solutions. We establish the general formalism for the stability 
considerations of these solutions and give examples of stability diagrams. Finally, we show that the 
effective Peierls-Nabarro barrier for the pulse-like soliton solution is zero. 

PACS numbers: 

I. INTRODUCTION. 

The discrete nonlinear Schrodinger equation (DNLSE) finds widespread use in physics due to its very general 
nonlinear character. It arises in the context of the propagation of electromagnetic waves in optical waveguides [1], 
and it also appears in the study of Bose-Einstein condensates in optical lattices [2|]. Recently we have studied the 
exact soliton solutions and their stability for the one-dimensional DNLSE with a saturable nonlinearity [3| . We were 
also able to obtain staggered and short period solutions of this equation Q as well as to generalize our results to 
arbitrarily higher-order nonlinearities Q|. 

Two-dimensional periodic lattices with a saturable nonlinearity in the Schrodinger equation have been experimen- 
tally realized in photorefractive materials Q . Solitons have also been observed in these crystals Q| ■ Localized 
traveling wave solutions that exist only for finite velocities have been computed in this case The question of 



10]. Here we show that it is also possible to 



discrete soliton mobility in these systems has been addressed as well 
obtain exact periodic and pulse-like soliton solutions for the two-dimensional DNLSE with a saturable nonlinearity. 
We then study the stability of these solutions as a function of the parameters of the equation. We also show that 
similar to the one dimensional case, the effective Peierls-Nabarro barrier (i.e. the discreteness barrier for soliton 
motion) for the pulse-like soliton solutions is zero. In addition, we find several short period solutions. 
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II. TWO-DIMENSIONAL DISCRETE NONLINEAR SCHRODINGER EQUATION WITH 

SATURABLE NONLINEARITY. 



The equation we consider is the following asymmetric, DNLSE with a saturable nonlinearity in two dimensions 

^ #n,m + [£(0 + n -l,m) + £(0n,m+l + 0n,m-l)] + -TTTT ^"i'-T = • (!) 

at ' l+|<Pn,m| 

|£ — £| is a measure of the spatial asymmetry and v a measure of the nonlinearity. This equation can be derived from 

the Hamiltonian 

n.m r 



n,m—l 



(\<fcn+l,m - 0n,m| 2 + £\<j> 



n,m+l V^n,m| 



-[2(C + + v]\<t>n,m\ 2 +i/ln(l + |0 



and the equation of motion being derived from 



as? 



(2) 



(3) 



considering <p n . m and i0* m as conjugate variables. There are two conserved quantities for the field equation, Eq. 
([T]), the Hamiltonian H and the power (norm) P defined by 

N.M 



n,m— 1 

Note that the system is invariant under simultaneous interchange of £ and £, and n and m. 



(4) 



III. EXACT SOLUTIONS TO THE TWO-DIMENSIONAL EQUATION. 

Exact stationary solutions can also be obtained in the case of this two-dimensional discrete, asymmetric saturable 
nonlinear Schrodinger equation ([T]). We are looking for stationary solutions using the ansatz 

(5) 

and obtain from Eq. (JTJ) , the following difference equation 

(1 + «n.Jw«n,m + (1 + «n, m )[C(«n+l,m + " n -i,m) + £(«n,m+l + «n,m-l)] + ^ u «,m = ■ ( 6 ) 

Following Ref. Q we immediately find two different types of solutions. One that is symmetric in n and m and one 
that only depends on n (and by symmetry one that only depends on m). 



Symmetric case: 
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If one chooses 

« = -!/, (7) 

then 

<, m = ^§jM/3(n + m + 6i),k), (8) 



is a solution if k is chosen to fulfill 



C + C cn 2 (/3,fc) ' ^ N p 

Similarly 

sn(/3,fc) 



i/ _ dnQ8,fc) 2K(fc) 



is a solution provided 



<m ^ fc dn( y fc) cn(/3(n + m + fr), fc) , (10) 



»--o c f: h l /3=^. (id 



C + e dn^(/3,fc) iV p 
Here fc is the elliptic modulus (the elliptic parameter m — k 2 [11]) of the Jacobi elliptic functions sn(x, &:), cn(a;, k), 
and dn(x, k) and K(k) is the complete elliptic integral of the first kind [lj], [ijj]. The integer N p denotes the spatial 
period of the system. N and M in Eq. @ must be chosen as multiples of N p . The two solutions have a common 
pulse- like limit for k — > 1 (and AT p — >• oo), 

u s n m = sinh(/?)sech[/3(n + m + Si)] , (12) 

which is a solution if (3 fulfills 

-^i- = 2cosh(/3). (13) 

By symmetry, a change of sign of n (or m) in Eqs. (|8|). (jlpp . and (1121) will give solutions with exactly the same 
properties. 

Asymmetric case: 

If one chooses for the frequency 

u = -v-2£, (14) 

then 



(15) 



will be a solution for k satisfying 



Similarly, we also have a cn solution. 



provided 



v dnQ8,fc) 2K(k) 

C =2 ^nW P= ^- (16) 



= k ^{^) Cn{ P {n + Sl) > k) > (17) 



1/ cn(/3,fc) 4K(k) , N 

c= 2 ^' *~sr- <18) 



Here N — N p and M can be any integer. 

Again both these solutions approach the pulse solution in the limit k — > 1 and N p — > oo 

u a n s m = sinh(/3)sech[/3(n + Si)] , (19) 

provided 

^ = 2cosh(^). (20) 

Another asymmetric solution appears if we interchange n and m and £ and £ (M = N p and N any integer). We 
note that, in all cases, changing the sign of £ is equivalent to staggering the solution in the m direction ((— l)™ 1 as 
an amplitude factor) and changing the sign of £ is equivalent to staggering the solution in the n direction [4(]. 

The described solutions are in some sense direct generalizations of our earlier results for the one-dimensional version 
of the saturable nonlinear Schrodinger equation, since they remain spatially uniform along a specific direction in 
space. Such solutions are well-known for nonlinear partial differential equations and are often, in such a continuum 
setting, referred to as line solutions or line solitons. However, in continuum settings, such solutions are generally 
not stable, because the extra dimension now allows for an entire set of new instability modes to come into play. 
In our discrete case, however, we shall demonstrate that these solutions indeed can be stable in certain cases. The 
fact that these solutions have infinite extension along one of their dimensions probably renders them less physically 
important. However, their stability analysis does, as we will demonstrate, give detailed insight into the intricate 
stability mechanisms of this nonlinear system. Specifically, are the parameters {v, £), which control the stability, 
directly related to materials properties such as the change in refractive index of the crystal. 

IV. STABILITY OF THE SOLUTIONS. 

In order to study the linear stability of these exact solutions u J n m (J is "s" (symmetric) or "as" (asymmetric)) we 
introduce the following expansion around the exact solution 

<i>n, m (t) = < m e-^ t+5 > + 5u n , m (t)e-^ t+ V (21) 



applied in a frame rotating with frequency u> of the solution. Substituting (f2"Tj) into the field equation, Eq. ([!}, and 
retaining only terms linear in the deviation, <5«.„ >m , we get 

^u„ jTO +C(ou n +i +ou J+4(du i+du i)+ w _| — a u +- — = 0. 

V (l + K,m|T / (1 + | 

(22) 

We continue by splitting the deviations bu n , m into real parts dun,ln and imaginary parts 8un]m {8u n ^ m — 5un}m + 
i5un,m) and introducing the two real vectors 

SIT = {SuQj = {SU^}, and SU* = {6u% m } = {51%% (23) 

where the pair of indices m,n are replaced by a single index J via: J = n + (to — 1)-/V p . By introducing the real 
matrices A = {Aj j*} and B = {Bj ji} defined by 



Aj,j> = Wj,J>+i + Sj.J'-i) + t(5j,j, +Np + 5j,j,- Np ) + [ u, + H< ;" J2 | 3 J +l ^ 2 m ' 2) ] Sj,j>, (24) 

(1 T | U n ,m | ) 



V (l + |Un,m| 2 )/ 

where J' ± 1 and J' ± AT p in the Kronecker 6 means: J' ± 1 mod N p and J' ± iV p mod _/V p to ensure periodic boundary 
conditions, Eq. (|2"2"j) becomes 

- + A6U r = 0, and <5t> + BSW = 0. (26) 

Combining these first order differential equations we get: 

5U* + AB8U i = 0, and 5U r + BASU r = 0. (27) 

The two matrices A and B are symmetric and have real elements. However, since they do not commute AB 
and BA = (AB) tr are not symmetric. AB and BA have the same eigenvalues, but different eigenvectors. The 
eigenvectors for each of the two matrices need not be orthogonal. The eigenvalue spectrum {7} of the matrix AB 
(or BA) determines the stability of the exact solutions. If {7} contains negative eigenvalues then the solution is 
unstable. The eigenvalue spectrum always contains two eigenvalues which are zero. These eigenvalues correspond 
to the translational invariance in space and time (represented by <5i and 5). The given solutions are unstable for 
most of the parameter space (C:£, ^)- From Ref. [3j we know that £ = generally leads to stable solutions. In 
determining the stability of the solutions it is useful to note that the rescaling transformation (£, £, v) —¥ a(£,£, z^) 
(oj — > acu) changes the eigenvalues by a 2 and therefore it does not affect the stability (i.e. the sign of eigenvalues) 
of the solutions. Therefore the three-dimensional parameter space (£, £, v) can be significantly reduced (into a 
two-dimensional parameter space) as far as stability considerations are concerned. The nonlinearity parameter 
v 7^ 0, separates the three-dimensional parameter space into two disconnected equivalent ones for v > and v < 0, 
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respectively. The rescaling transformation with a — — 1 interchanges the two equivalent half spaces. So we need only 
to consider positive v. From here on we treat the symmetric and the asymmetric cases separately. 
Stability of the symmetric case: 

Since v > 0, Eqs. ©, (fTTT) . and (| 13[) require ( + £ > and therefore ( > — £. Further, we can always choose 
C > £ due to the inter-changeability of C and £. We therefore have — £ < £ < £ or applying the scaling condition 
— 1 < £/£ < 1. This means that the stability of the entire parameter space can be mapped out onto the much 
smaller parameter space (1,£, v), where — 1 < £ < 1 (and v > 0). In Fig. Q]we illustrate the stability analysis 




-2 -1.5 -1 -0.5 

\ 

FIG. 1: Stability analysis for the dn solution, Eq. {[Sj . Lowest eigenvalues for \u\ = 1.5 (thick line, black online), \v\ = 2 
(medium thick line, blue online), and \v\ = 3 (thin line, red online). Dashed (solid) lines indicate negative (positive) values of 
v. Stability occurs when the lowest eigenvalue is zero. The entire existence interval is shown for positive v. Negative v results 
can be obtained by rescaling the results for positive v. The remaining parameters are: £ = 1 and N p = 8. 

for the dn solution, Eq. ©, by showing the lowest eigenvalue as a function of £ for several values of v. Stability 
occurs whenever the lowest eigenvalue is zero. The entire existence regime is illustrated for positive v and a few 
windows of stability can be seen. It is important to note that the results for negative values of v (£ < —1) are 
superfluous as they can be obtained by rescaling the results for positive v. To demonstrate this we have for £ = 1: 
(1, — — ( — 1/ 1^| , 1, M/|£|) — > (1, — 1/|£|) where the last step follows by the inter-changeability of the 

two-coupling parameters. This shows that the v < regime can be mapped onto the v > regime. Generally, 
stability is only observed for £ < 0. This means that stability for the symmetric dn solutions can only be achieved 
when £ and C have opposite signs. Recalling the equivalence between staggered solutions and sign changes of £ or 
C, another way to view this result is that the symmetric dn solutions must be staggered in one dimension when £ 
and £ are both positive in order to be stable. This is clearly a property arising from the discreteness that cannot be 
achieved in a continuum system. 
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FIG. 2: General stability diagram for symmetric dn (left) and cn (right) solutions of Eq. (JSJ and Eq. (|10p . respectively. In 
the grey region no solution of the given kind exists. In the white regions solutions are unstable, while the black area indicates 
stable solutions. The vertical line £ = 0, which represents the one-dimensional stability result, is part of the black area. The 
rough or jagged appearance of parts of the stability boundary is due to limited numerical resolution rather than an intrinsic 
feature of the problem. Parameters are: £ = 1, and N p = 8. 

Assembling results like those shown in Fig. Q] for a range of v and £ values we arrive at Fig. [2] where stability 
diagrams for both the dn and the cn solutions to Eq. ((5J and Eq. (|T0|) . respectively, are shown. The grey region 
indicates that the solutions do not exist, whereas the black (white) regions indicate the existence of stable (unstable) 
solutions. These two stability diagrams have a very similar structure. However, the cn solutions are always stable 
in the proximity of the existence boundary marked by the grey area. This property is related to the fact that the 
amplitude (which is oc k, see Eq.(fT0])) of this solution vanishes at the existence boundary where k — 0. Also, we note 
that the common pulse solution corresponds to the corner close to (£, v) = (—1, 0). 

We have looked at other values of N p and find that for larger N p , the stability diagram has similar features. The 
pulse- like solution only exists in the limit N p — ¥ oo, and here it has the same stability properties as the dn and cn 
solutions for k — >■ 1. Therefore, our analysis indicates that the symmetric pulse- like solutions are stable for small v 
and £ ~ — 1. 

Stability of the asymmetric case: 

We proceed almost as in the symmetric case. We still have v > 0, therefore from Eqs. (fl6]l. (p~8]> . and (|20j) we 
have £ > 0. So here the three-dimensional parameter space can be reduced to (1,£, v), where — oo < £ < oo (and 
v > 0). Illustrations of the stability diagrams are given for the asymmetric dn and cn solutions, Eq. (|15[) and Eq. 
(|17p. respectively in Fig. [3] for £ = 1 and N p = 8. Again the two stability diagrams have a very similar structure 
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FIG. 3: General stability diagram for the asymmetric dn (left) and cn (right) solutions of Eq. (|15[) and Eq. JTJJ, respectively. 
In the grey region no solution of the given kind exists. In the white region solutions are unstable, while the black area 
indicates stable solutions. The vertical line £ = 0, which represents the one-dimensional stability result, is part of the black 
area. Parameters are: £ — 1, and N p = 8. 

except that the parameter space for stability of the asymmetric solution is much larger. As in the symmetric case, 
even in the asymmetric case, the cn solutions are always stable in the proximity of the existence boundary marked 
by the grey area. Here we note that the common pulse solution corresponds to large v and £. 



V. PEIERLS-NABARRO BARRIER FOR THE PULSE SOLUTION 



We would now like to show the absence of Peierls-Nabarro barrier for the pulse solution. However, we must 
remember that since both power P and Hamiltonian H are constants of motion, one must compute the energy 
difference between the solutions when Si — and Si = 1/2 in such a way that the power P is same in both the cases. 

For the pulse solution obtained above, the power P is given by 



p = Yj l<^™| 2 = s inh 2 (^) Yj sech2 [P( n + m + S i)} ■ 

n,m — — oo n,m — — oo 

This double sum can be evaluated using the single sum result 

Y 2 - L 2, 0/ r x, 2 2K(k)E(k) (2K{k)\ 2 , 2 . . , 
Y sech 2 [/3(n + 5 2 )] = - + ^ + ( — ) dn 2 [2S 2 K(k),k 



71— — OO 



the above P is given by 

oo 

P= ]T sinh 2 (/3) 



m— — oo 



2 2K(k)E(k) f2K(k)\' 



dn 2 [2(m + 6i)K(k),k] 



(28) 



(29) 



(30) 
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Note, only the last term on the rhs is to dependent. In these equations E(k) is the complete elliptic integral of the 
second kind. 

Let us now discuss the computation of the Hamiltonian H . Clearly, for the pulse solution obtained above, H as 
given by Eq. @ takes the form 



*= E 



n,m=— go 



- vP + vln[l + sinh 2 (/3)sech 2 (/3[n + m + 5i])] 



-2sinh 2 (/3)[sech(/3[n + 1 + to + <5i])sech(/3[n + m + 5 1 })(( + £)] 



Again we use the single sum results to evaluate the double sum, i.e. 



J2 [sech[£(n + 1 + <5 2 )]sech[/3(n + S 2 )}] = 



n— — oc 



sinh(/3) ' 



(31) 



(32) 



J2 H 1 + sinh 2 (/3)sech 2 (^[n + 5 2 ])} = 2/3 , 



n— — oo 



the above H is given by 



vP + 2v/3 - 4(C + sinh(/3) 



(33) 



(34) 



We thus note that for a given power P (which contains a sum over to), H is indeed independent of 5\, i.e. the 
Peierls-Nabarro barrier is indeed zero for the pulse solution. The same holds true for the asymmetric solution. 



VI. SHORT PERIOD SOLUTIONS 

Recently we obtained short period solutions to the one-dimensional saturable DNLSE Q . These short period ( N) 
solutions, in the one-dimensional case, can be written in the following compact form (coming from equally distributed 
points on a circle so that projection on the x-axis should only be or ±a): 

/2?m 



MAr(n) 



COS((y3Af) \ N 

where (pi =(p 2 =(fi4 S = 0, (4s is the stable period 4) and ip4 = tps = ipe = -| . 



u N {n + 1) + u N (n - 1) = 



cos(^jv) 



2vr(n + l) \ /27r(n-l) 
cos | — h (pjsr + cos — h ifN 



N 



N 



2a ( 2im \ 2tt , . . 2vr 

__ cos ^_ + ^ j cos _ = 2M n) cos -. 



For lo we get 



oj = — 2<^cos 



2tt va 1 
~~N ~ 1 + a 2 



(35) 

(36) 
(37) 

(38) 
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Assuming that in the two-dimensional case, the solution is a product of the two one-dimensional solutions (with 
period N and period M), i.e. 

i s a ( \ / 27rm \ . . 

U N , M (n, m) = r r COS — - -I- tp N COS — — + tfM , (39) 

cos(ip N )cos(ip M ) \ N J \ M J 



we get: 



2n _ 2n va 2 , . 

uj = -2Ccos— - 2^ cos— - -. (40) 



N M 1 + a 



VII. CONCLUSIONS 



We have given analytical expressions for the solutions to the two-dimensional discrete nonlinear Schrodinger 
equation with saturable nonlinearity which arises in photorefractive crystals Q-l^. Due to their infinite extension 
along one of their dimensions, these solutions are not very physically meaningful but it is very rare that solutions to 
discrete nonlinear two-dimensional problems can be described in closed form using standard mathematical functions 
as we have done here. This feature of the solutions is physically significant because it provides an opportunity for 
in-depth scrutiny and understanding that is not usually available in a nonlinear physical system. These solutions are 
closely related to the previously derived 3] solutions to the corresponding one-dimensional equation. However, in 
contrast to what one may expect based on intuition derived from similar nonlinear partial differential equations, we 
have shown that these solutions are linearly stable in certain regions of the parameter space. Specifically, we have 
observed that the asymmetric versions of these solutions lead to a very intricate stability diagram. We have shown 
that the symmetric dn solution is stable in certain regions of the parameter space provided it is staggered in one 
dimension. However, the symmetric cn solution as well as the asymmetric cn and dn solutions are stable in certain 
regions of parameter space both when they are non-staggered or if they are staggered in one dimension. The finding 
that nonlinear waveforms in two-dimensional photorefractive materials best achieve stability in the presence of phase 
asymmetry between the two spatial directions is crucial because the photonic lattices that represent the physical 
realization of Eq.([T|) tend to naturally possess this property 14|. Finally, we found that the Peierls-Nabarro barrier 



for the pulse solution is zero. An understanding of the mobility of these exact discrete two-dimensional solutions 
remains an important issue [lo| . 
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